function gesol = PriceConsClaimUnlevered(parms,gesol_in,ERRTOL)
% P = C + E[M'P'], cum dividend price

    MAXITERS = 10000;

    gesol = gesol_in;
    
    NN = parms.NN; Nx = parms.Nx; Nz = parms.Nz;
    
    C_Nxz = gesol.C(:);
    P_Nxz_old = reshape(C_Nxz,[NN Nx Nz]);
    P_Nxz = zeros(NN,Nx,Nz);
    
    % laws of motion
    x_next_Nxzz = gen_law_of_motion_x(parms);
    Nnext_Nxz = gen_law_of_motion_N(parms,gesol);
    
    [N_meshgrid,x_meshgrid] = meshgrid(parms.grid_N(:),parms.grid_x(:)); % needed for interp2
    P_next = zeros(NN*Nx*Nz,Nz); % P'=P'((N,L,z),z')
    
    StateTransitionProbs_NLzz = reshape(permute(repmat(parms.StateTransitionProbs,[1 1 NN Nx]),[3 4 1 2]),...
        [NN*Nx*Nz, Nz]);
    ADPrice_NLzz = StateTransitionProbs_NLzz.*reshape(gesol.spd1period,[NN*Nx*Nz, Nz]);
    
    bContinue = true;
    iter = 0;
    
    while bContinue
        
        for iz_next = 1:parms.Nz
            P_next(:,iz_next) = interp2(N_meshgrid,x_meshgrid,P_Nxz_old(:,:,iz_next)',...
                Nnext_Nxz(:),x_next_Nxzz(:,iz_next));
        end
        
        P_Nxz(:) = C_Nxz + sum(ADPrice_NLzz.*P_next,2);
        
        VFIerr = max(abs(P_Nxz(:) - P_Nxz_old(:)));
        P_Nxz_old(:) = P_Nxz;
        
        iter = iter + 1;
        
        if VFIerr<=ERRTOL || iter>=MAXITERS
            bContinue = false;
        end
        
    end
    
    gesol.stock.price_cons_unlev = reshape(P_Nxz,[NN Nx Nz]);
    gesol.stock.note = 'P = C + E[M''*P''] is the cum-dividend price.';
    gesol.stock.VFIerr = VFIerr;

end